home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / DBio / DSequence.cpp < prev    next >
Encoding:
Text File  |  1995-12-17  |  34.3 KB  |  1,488 lines  |  [TEXT/R*ch]
  1. // DSequence.cp 
  2. // d.g.gilbert
  3.  
  4.  
  5.  
  6. #include "DSequence.h"
  7. #include "DSeqFile.h"
  8. #include "DREnzyme.h"
  9. #include <ncbi.h>
  10. #include <dgg.h>
  11. #include <DUtil.h>
  12. #include <DFile.h>
  13. #include "ureadseq.h"
  14.  
  15.  
  16. char* DSequence::kConsensus = "Consensus";
  17. char  DSequence::indelHard = '-';
  18. char  DSequence::indelSoft = '~';
  19. char  DSequence::indelEdge = '.';
  20.  
  21.  
  22. Local char* gAminos        =    "ABCDEFGHIKLMNPQRSTVWXYZ*_.-~?";
  23. Local char* gIubbase    = "ACGTUMRWSYKVHDBXN_.-*~?";     //did include ".", for what?
  24. Local char* gPrimenuc    = "ACGTU_.-*~?";
  25. Local char* protonly    = "EFIPQZ";
  26. Local char* stdsymbols= "_.-*?";
  27. Local char* allsymbols= "_.-*?!=/**/[]()!&#$%^&=+;:/|`~'\"\\";
  28. Local char* seqsymbols= allsymbols;
  29. Local Boolean kForWriting    = TRUE;
  30. Local Boolean    kForReading    = FALSE;
  31.  
  32. enum { kSearchNotFound = -1, kMaskOK = 'K', kMaskEmpty = -1 };
  33.  
  34.  
  35.     
  36. long DSequence::NucleicBits(char nuc)
  37. {
  38.     //char nuc= toupper(nuc);
  39.     switch (nuc) {
  40.         case 'a': case 'A': return kMaskA;  
  41.         case 'c': case 'C': return kMaskC;   
  42.         case 'g': case 'G': return kMaskG;   
  43.         case 't': case 'T': return kMaskT;   
  44.         case 'u': case 'U': return kMaskU;  // ==  ( kMaskT | kMaskExtra); 
  45.         case 'y': case 'Y': return (kMaskC | kMaskT);  
  46.         case 'k': case 'K': return (kMaskG | kMaskT);  
  47.         case 's': case 'S': return (kMaskG | kMaskC);  
  48.         case 'r': case 'R': return (kMaskG | kMaskA);  
  49.         case 'm': case 'M': return (kMaskA | kMaskC);  
  50.         case 'w': case 'W': return (kMaskA | kMaskT);  
  51.         case 'h': case 'H': return (kMaskA | kMaskC | kMaskT);  
  52.         case 'b': case 'B': return (kMaskG | kMaskC | kMaskT);  
  53.         case 'v': case 'V': return (kMaskG | kMaskC | kMaskA);  
  54.         case 'd': case 'D': return (kMaskG | kMaskT | kMaskA);  
  55.         // case '-': return kMaskIndel;   
  56.         case ' ': case '_': return 0;  //? spacers
  57.         case  0 : return 0;  
  58.         // case '.': //indelEdge     
  59.         case 'n': case 'N': return kMaskNucs; 
  60.         default : 
  61.             {
  62.             if (nuc == indelHard || nuc == indelSoft)  return kMaskIndel;
  63.             else if (nuc == indelEdge) return kMaskNucs;
  64.             else return 0; //kMaskNucs;  //?? match all or match none ??
  65.             }
  66.         }
  67. }
  68.  
  69.  
  70. Boolean DSequence::IsNucleicMatch( long xNucBits, long yNucBits)
  71. {
  72.         // note: Indel/spaces match nothing, including other indel/space ...
  73.     return 0 != ( (kMaskNucs & xNucBits) & (kMaskNucs & yNucBits));
  74. }
  75.  
  76.  
  77. char DSequence::NucleicConsensus( long xNucBits, long yNucBits)
  78. { 
  79.     char nuc; 
  80.     Boolean isrna, indel;
  81.  
  82.     if (xNucBits==0 || yNucBits==0) 
  83.         return ' ';  //?? either a space yields space ?
  84.     else {
  85.         indel= (kMaskIndel == xNucBits || kMaskIndel == yNucBits);
  86.         isrna= (kMaskU == xNucBits || kMaskU == yNucBits);
  87.         switch ( (kMaskNucs & xNucBits) | (kMaskNucs & yNucBits)) {
  88.             case kMaskA: nuc= 'A'; break;
  89.             case kMaskC: nuc= 'C'; break;
  90.             case kMaskG: nuc= 'G'; break;
  91.             case kMaskT: 
  92.                 if (isrna) nuc= 'U'; else nuc= 'T';
  93.                 break;
  94.                     //? do other ambig codes ?
  95.             default:    nuc= '.'; break; //? use this as Consensus other ?
  96.             }
  97.         if (indel) nuc= tolower(nuc); //?
  98.         }
  99.     return nuc;
  100. }
  101.  
  102.  
  103.  
  104. void DSequence::NucleicComplement( Boolean isrna, char* fromSeq, char* toSeq, long len)
  105. {
  106. // note: fromSeq may == toSeq 
  107.     register char b, c;
  108.     Boolean islow;
  109.      while (len--) {
  110.         b= *fromSeq++;
  111.         islow= islower(b);
  112.         if (islow) b= toupper(b);
  113.         switch (b) {
  114.             case 'C': c= 'G'; break;
  115.             case 'G': c= 'C'; break;
  116.             case 'A': if (isrna) c= 'U'; else c= 'T'; break;
  117.             case 'T':
  118.             case 'U': c= 'A'; break;
  119.             case 'R': c= 'Y'; break;
  120.             case 'Y': c= 'R'; break;
  121.             case 'M': c= 'K'; break;
  122.             case 'K': c= 'M'; break;
  123.             case 'S': c= 'S'; break;  // was mistaken 'w'
  124.             case 'W': c= 'W'; break;  // was mistaken 's'
  125.             case 'H': c= 'D'; break;
  126.             case 'D': c= 'H'; break;
  127.             case 'B': c= 'V'; break;
  128.             case 'V': c= 'B'; break;
  129.             default : c= b;   break;
  130.             }
  131.         if (islow) c= tolower(c);
  132.         *toSeq++= c;
  133.         }
  134. }
  135.  
  136.  
  137.  
  138. const char* DSequence::Amino123( char amino1)
  139. {
  140.  switch (toupper( amino1)) {
  141.      case 'A': return "Ala";
  142.      case 'M': return "Met";
  143.      case 'R': return "Arg";
  144.      case 'F': return "Phe";
  145.      case 'N': return "Asn";
  146.      case 'P': return "Pro";
  147.      case 'D': return "Asp";
  148.      case 'S': return "Ser";
  149.      case 'C': return "Cys";
  150.      case 'T': return "Thr";
  151.      case 'Q': return "Gln";
  152.      case 'W': return "Trp";
  153.      case 'E': return "Glu";
  154.      case 'Y': return "Tyr";
  155.      case 'G': return "Gly";
  156.      case 'V': return "Val";
  157.      case 'H': return "His";
  158.      case 'I': return "Ile";
  159.      case 'B': return "Asx";
  160.      case 'L': return "Leu";
  161.      case 'Z': return "Glx";
  162.      case 'K': return "Lys";
  163.      case 'X': return "Xaa";
  164.      case '*': return "End";
  165.      case ' ': return "   ";
  166.      default : return "???"; //??
  167.      }
  168. }
  169.  
  170.  
  171. char DSequence::Amino321( char* amino)
  172. {
  173.     char    aa[4];
  174.     aa[0]= toupper(amino[0]);
  175.     aa[1]= tolower(amino[1]);
  176.     aa[2]= tolower(amino[2]);
  177.     aa[3]= 0;
  178.     
  179.              if (!StrCmp(aa,"Ala")) return 'A';
  180.     else if (!StrCmp(aa,"Met")) return 'M';
  181.     else if (!StrCmp(aa,"Arg")) return 'R';         
  182.     else if (!StrCmp(aa,"Phe")) return 'F';
  183.     else if (!StrCmp(aa,"Asn")) return 'N';       
  184.     else if (!StrCmp(aa,"Pro")) return 'P';
  185.     else if (!StrCmp(aa,"Asp")) return 'D';    
  186.     else if (!StrCmp(aa,"Ser")) return 'S';
  187.     else if (!StrCmp(aa,"Cys")) return 'C';         
  188.     else if (!StrCmp(aa,"Thr")) return 'T';
  189.     else if (!StrCmp(aa,"Gln")) return 'Q';       
  190.     else if (!StrCmp(aa,"Trp")) return 'W';
  191.     else if (!StrCmp(aa,"Glu")) return 'E';    
  192.     else if (!StrCmp(aa,"Tyr")) return 'Y';
  193.     else if (!StrCmp(aa,"Gly")) return 'G';          
  194.     else if (!StrCmp(aa,"Val")) return 'V';
  195.     else if (!StrCmp(aa,"His")) return 'H';
  196.     else if (!StrCmp(aa,"Ile")) return 'I';       
  197.     else if (!StrCmp(aa,"Asx")) return 'B';
  198.     else if (!StrCmp(aa,"Leu")) return 'L';          
  199.     else if (!StrCmp(aa,"Glx")) return 'Z';
  200.     else if (!StrCmp(aa,"Lys")) return 'K';           
  201.     else if (!StrCmp(aa,"End")) return '*';           
  202.     else if (!StrCmp(aa,"   ")) return ' '; //??           
  203.     else return 'X';
  204. }
  205.  
  206.  
  207. char* DSequence::Nucleic23Protein( char* nucs, long nbases)
  208. //!! Need to do this w/ NucBits so ambiguity codes work !!
  209. //!! The CodonTable.codon == acodon should be IsNucleicMatch( bits(codon), bits(acodon)) 
  210. {
  211.     long    i;  //borland 16bit don't like longs for iterators !?
  212.     char  *ap, *kp, *aminos;
  213.  
  214.     if (DCodons::NotAvailable()) return NULL; // (gCodonTable[0].amino == 0)
  215.  
  216.     aminos= StrDup( nucs);
  217.     for (i= 0, ap= aminos; i<nbases; i++) {
  218.         register char c= toupper( *ap);
  219.         if (c=='U') c= 'T';
  220.         *ap++= c;
  221.         }
  222.  
  223.     for (i= 0, ap= aminos, kp= aminos; i<nbases; i++) {
  224.         for (long j=0; j<64; j++)
  225.             if ( MemCmp(DCodons::codon(j), ap, 3) == 0) { //gCodonTable[j].codon
  226.                 *kp++= DCodons::amino(j); // gCodonTable[j].amino;
  227.                 break;
  228.                 }
  229.         }
  230.     *kp= '\0';
  231.     aminos= (char*) MemMore( aminos, kp - aminos + 1);
  232.     return aminos;
  233. } 
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243. // DSequence -----------------
  244.  
  245.  
  246. DSequence::DSequence()
  247. {
  248.     fBases= (char*) MemGet(1,true); 
  249.   fMasks= NULL; //(char*) MemGet(1,true); 
  250.   fMasksOk= false;
  251.     fInfo    = (char*) MemGet(1,true);     
  252.     fName = StrDup("Untitled");         //!? need unique name, eg. Noname-1,noname-2...
  253.     fKind = kOtherSeq;         //? user default
  254.     fLength        = 0;
  255.     fChecksum    = 0;
  256.     fIndex        = 0; //dummy
  257.     fModTime    = gUtil->time();
  258.     fChanged    = TRUE;
  259.     fDeleted    = FALSE; 
  260.     fOpen            = FALSE;
  261.     fSearchRec.targ= NULL;
  262.     fOrigin    = 0;
  263.     fSelStart= 0;
  264.     fSelBases= 0;
  265. }
  266.  
  267.  
  268. DSequence::~DSequence()  
  269. {
  270.     MemFree(fBases);
  271.     MemFree(fMasks);
  272.     MemFree(fInfo);
  273.     MemFree(fName);
  274.     MemFree(fSearchRec.targ);
  275. }
  276.  
  277.  
  278. DObject* DSequence::Clone() // override 
  279. {    
  280.     DSequence* aSeq= (DSequence*) DObject::Clone();
  281.     aSeq->fBases= StrDup( fBases);
  282.     aSeq->fMasks= NULL; aSeq->SetMasks( fMasks); //aSeq->fMasks= StrDup( fMasks);
  283.     aSeq->fInfo=  StrDup( fInfo);
  284.     aSeq->fName=  StrDup( fName);
  285.     aSeq->fSearchRec.targ= StrDup( fSearchRec.targ);
  286.     return aSeq;
  287. }
  288.  
  289. void DSequence::CopyContents( DSequence* fromSeq)
  290. {
  291.     SetBases( fromSeq->fBases); 
  292.     SetMasks( fromSeq->fMasks); 
  293.     SetInfo( fromSeq->fInfo); 
  294.     SetName( fromSeq->fName); 
  295.     fKind            = fromSeq->fKind;          
  296.     fLength        = fromSeq->fLength;
  297.     fChecksum    = fromSeq->fChecksum;
  298.     fModTime    = fromSeq->fModTime;
  299.     fChanged    = fromSeq->fChanged;
  300.     fDeleted    = fromSeq->fDeleted; 
  301.     fOpen            = fromSeq->fOpen;
  302.     fIndex        = fromSeq->fIndex;
  303.     fSearchRec= fromSeq->fSearchRec;
  304.     fSearchRec.targ= StrDup( fromSeq->fSearchRec.targ);
  305. }
  306.  
  307. void DSequence::SetBases( char* theBases)
  308. {
  309.     if (fBases) MemFree(fBases);
  310.     fBases = StrDup(theBases);  
  311.     fLength= StrLen(fBases);
  312.     fMasksOk= false; // assume SetBases allways called before SetMask...
  313. }
  314.  
  315. void DSequence::UpdateBases( char* theBases, long theLength)
  316. {
  317.     fBases =  theBases;  
  318.     fLength= theLength;
  319.     fMasksOk= false; // assume this allways called before SetMask...
  320. }
  321.  
  322.  
  323. void DSequence::SetMasks( char* theMasks)
  324. {
  325.     if (fMasks) MemFree(fMasks);
  326.     if (theMasks) {
  327.         fMasks = StrDup(theMasks);  
  328.         FixMasks();
  329.         }
  330.     else {
  331.         fMasks= NULL;
  332.         fMasksOk= false;
  333.         }
  334. }
  335.  
  336. void DSequence::FixMasks()
  337. {
  338.     ulong i, mlen= 0;
  339.     if (!fMasks) { 
  340.         mlen= 1+fLength;
  341.         fMasks= (char*) MemGet(mlen,true);  
  342.         for (i=0; i<mlen; i++) fMasks[i] |= 0x80; // set high bit so Strxxx sees it as non-null
  343.         fMasks[mlen-1]= 0;
  344.         }
  345.     else { 
  346.         mlen= 1+StrLen(fMasks);
  347.         if (fLength+1 != mlen) {
  348.             fMasks= (char*) Nlm_MemExtend( fMasks, fLength+1, mlen);
  349.             mlen= fLength+1; 
  350.             for (i=0; i<mlen; i++) fMasks[i] |= 0x80; 
  351.             fMasks[mlen-1]= 0;
  352.             } 
  353.         }
  354.     fMasksOk= true;
  355. }
  356.  
  357. short DSequence::MaskAt(long baseindex, short masklevel)
  358. {
  359.     short b = kMaskEmpty;
  360.     if (fMasks && fMasksOk && baseindex>=0 && baseindex<fLength) {
  361.         b= (unsigned char) fMasks[baseindex];
  362.         switch (masklevel) {
  363.             case 0: b &= 0x7f; break; // return full mask - top (0x80) bit
  364.             case 1: b &= 1;  break;  
  365.             case 2: b &= 2;  break;
  366.             case 3: b &= 4;  break;
  367.             case 4: b &= 8;  break;
  368.             default: b = kMaskEmpty; break;
  369.             }
  370.         }
  371.     return b;
  372. }
  373.  
  374.  
  375. inline void SetMaskByte( char& b, short masklevel, short maskval)
  376. {
  377.     switch (masklevel) {
  378.         case 1: if (maskval) b |= 1; else b &= ~1; break;
  379.         case 2: if (maskval) b |= 2; else b &= ~2; break;
  380.         case 3: if (maskval) b |= 4; else b &= ~4; break;
  381.         case 4: if (maskval) b |= 8; else b &= ~8; break;
  382.         // reserve bits 5, 6, 7 & 8 for now 
  383.         // -- 7 & 8 may go unused to store data in printchar form
  384.         default: break;
  385.         }
  386.     //if (b) masklevel--; // debugging b val
  387. }
  388.  
  389. void DSequence::SetMaskAt(long baseindex, short masklevel, short maskval)
  390. {
  391.     if (fMasks && fMasksOk && baseindex>=0 && baseindex<fLength)  
  392.         ::SetMaskByte(fMasks[baseindex],masklevel, maskval);
  393. }
  394.  
  395.  
  396. void DSequence::SetMask(short masklevel, short maskval)
  397. {
  398.     long  baseindex;
  399.     if (fMasks && fMasksOk) 
  400.      for (baseindex=0; baseindex<fLength; baseindex++) 
  401.         ::SetMaskByte(fMasks[baseindex], masklevel, maskval);
  402. }
  403.  
  404.  
  405. inline void FlipMaskByte( char& b, short masklevel)
  406. {
  407.     switch (masklevel) {
  408.         case 1: b ^= 1;  break;  
  409.         case 2: b ^= 2;  break;
  410.         case 3: b ^= 4;  break;
  411.         case 4: b ^= 8;  break;
  412.         default:  break;
  413.         }
  414.     //if (b) masklevel--; // debugging b val
  415. }
  416.  
  417. void DSequence::FlipMaskAt(long baseindex, short masklevel)
  418. {
  419.     if (fMasks && fMasksOk && baseindex>=0 && baseindex<fLength)  
  420.         ::FlipMaskByte(fMasks[baseindex],masklevel);
  421. }
  422.  
  423. void DSequence::FlipMask(short masklevel)
  424. {
  425.     long  baseindex;
  426.     if (fMasks && fMasksOk) 
  427.      for (baseindex=0; baseindex<fLength; baseindex++)  
  428.         FlipMaskByte(fMasks[baseindex],masklevel);
  429. }
  430.  
  431.  
  432. void DSequence::SetInfo( char* theInfo)
  433. {
  434.     if (fInfo) MemFree(fInfo);
  435.     fInfo = StrDup(theInfo);  
  436. }
  437.  
  438. void DSequence::SetName( char* theName)
  439. {
  440.     if (fName) MemFree(fName);
  441.     fName = StrDup(theName);  
  442. }
  443.  
  444. void DSequence::SetKind( short theKind)
  445. {
  446.     fKind= theKind;  
  447. }
  448.  
  449. void DSequence::SetIndex( short theIndex)
  450. {
  451.     fIndex= theIndex;  
  452. }
  453.  
  454. char* DSequence::KindStr() const 
  455. {
  456.     switch (fKind) {
  457.         case kDNA             : return "DNA";
  458.         case kRNA             : return "RNA";
  459.         case kNucleic   : return "Nucleic";
  460.         case kAmino         : return "Protein";
  461.         default                    :
  462.         case kOtherSeq  : return "Unknown";
  463.         }
  464. }
  465.  
  466.  
  467. void DSequence::SetTime( unsigned long theTime)
  468. {
  469.     fModTime= theTime;
  470. }
  471.  
  472. void DSequence::SetOrigin( long theOrigin)
  473. {
  474.     fOrigin= theOrigin;
  475. }
  476.  
  477. void DSequence::Modified()  
  478. {
  479.     fModTime    = gUtil->time();
  480.     fChanged    = TRUE;
  481. }
  482.  
  483. void DSequence::Reformat(short format)
  484. {
  485.     // later...? 
  486. }
  487.  
  488. void DSequence::UpdateFlds()
  489. {
  490.     long         seqlen;
  491.     short     kind;
  492.     unsigned long check;
  493.     fLength= StrLen(fBases);
  494.     GetSeqStats(fBases, fLength, &kind, &check, &seqlen);
  495.     //fValidLength= seqlen; //?? is fLength == StrLen(fBases) or is it <= that, only valid bases??
  496.     fKind            = kind;
  497.     fChecksum    = check;
  498.     fChanged    = FALSE;
  499. }
  500.  
  501.  
  502.  
  503.  
  504.  
  505. void DSequence::SetSelection( long start, long nbases)
  506. {
  507.     fSelStart= start;
  508.     fSelBases= nbases;
  509. }
  510.  
  511. void DSequence::ClearSelection()
  512. {
  513.     fSelStart= 0;
  514.     fSelBases= 0;
  515. }
  516.  
  517. void    DSequence::GetSelection( long& start, long& nbases)
  518. {
  519.     if (fSelBases==0) {
  520.         start= 0;
  521.         nbases= fLength;
  522.         }    
  523.     else if (fLength>0) {
  524.         start= fSelStart;
  525.         nbases= Min( fSelBases, fLength - fSelStart);
  526.         }    
  527.     else {
  528.         start= 0;
  529.         nbases= 0;
  530.         }
  531. }
  532.  
  533.  
  534. DSequence* DSequence::CopyRange()
  535. {  
  536.             //note: start == 0 == 1st position 
  537.     DSequence* aSeq= (DSequence*) Clone();
  538.     if (fSelStart > 0 || fSelBases > 0) {
  539.         char* h= aSeq->fBases;
  540.         char* m= aSeq->fMasks;
  541.         Boolean domask= (fMasks && fMasksOk);
  542.  
  543.         long len= aSeq->fLength;
  544.         if (fSelStart > 0) {
  545.             len -= fSelStart;
  546.             MemMove( h, h + fSelStart, len);
  547.             if (domask) MemMove( m, m + fSelStart, len);
  548.             }
  549.         if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  550.         h= (char*) MemMore( h, len+1);
  551.         h[len]= '\0'; //((char*)h+len) = '\0';
  552.         aSeq->fBases= h;
  553.         if (domask) {
  554.             m= (char*) MemMore( m, len+1);
  555.             m[len]= '\0';    
  556.             aSeq->fMasks= m;
  557.             }
  558.         }
  559.     aSeq->UpdateFlds();
  560.     aSeq->Modified();
  561.     return aSeq;
  562. }
  563.  
  564.  
  565. DSequence* DSequence::CutRange()
  566. {
  567. /* aSeq.Cutrange means:  
  568.         (a) selected range (fSelStart..fSelStart+fSelBases) is copied into Result; &&
  569.         (b) selected range is removed from this
  570. */
  571.  
  572.     DSequence* aSeq= CopyRange();
  573.     if (fSelStart > 0 || fSelBases > 0) {
  574.         if (fSelStart<0) fSelStart= 0;
  575.         char* h= this->fBases;
  576.         char* m= this->fMasks;
  577.         Boolean domask= (fMasks && fMasksOk);
  578.         long len= this->fLength;
  579.         long top= len - fSelStart - fSelBases;
  580.         if (fSelBases > 0 && fSelBases < len) {
  581.             MemMove( h+fSelStart, h+fSelStart+fSelBases, top);
  582.             if (domask) MemMove( m+fSelStart, m+fSelStart+fSelBases, top);
  583.             }
  584.         len= fSelStart + top;
  585.         h= (char*) MemMore( h, len+1);
  586.         *((char*)h+len)= '\0';       //h[len]= '\0';
  587.         this->fBases= h;
  588.         if (domask) {
  589.             m= (char*) MemMore( m, len+1);
  590.             *((char*)m+len)= '\0';    
  591.             this->fMasks= m;
  592.             }
  593.         }
  594.     UpdateFlds();    
  595.     Modified();
  596.     return aSeq;
  597. }
  598.  
  599. void DSequence::InsertSpacers( long insertPoint, long insertCount, char spacer)
  600. { 
  601.     long i;
  602.     char *cp;
  603.     char* h= this->fBases;
  604.     char* m= this->fMasks;
  605.     Boolean domask= (fMasks && fMasksOk);
  606.     long len= this->fLength;
  607.     if (insertPoint<0) insertPoint= 0;
  608.     else if (insertPoint>len-1) insertPoint= len-1;
  609.     h= (char*) MemMore( h, len+insertCount+1);
  610.     cp = h+insertPoint;
  611.     MemMove( cp+insertCount, cp, len-insertPoint);
  612.     for ( i= 0; i<insertCount; i++)  *cp++= spacer;
  613.     h[len+insertCount]= '\0';
  614.     this->fBases= h;
  615.     if (domask) {
  616.         m= (char*) MemMore( m, len+insertCount+1);
  617.         cp = m+insertPoint;
  618.         char cmask= *cp;
  619.         MemMove( cp+insertCount, cp, len-insertPoint);
  620.         for ( i= 0; i<insertCount; i++)  *cp++= cmask;
  621.         m[len+insertCount]= '\0';
  622.         this->fMasks= m;
  623.         }
  624.     this->UpdateFlds();    
  625.     this->Modified();
  626. }
  627.  
  628.  
  629. DSequence* DSequence::PasteBases( long insertPoint, char* fromBases, char* fromMasks)
  630. { 
  631.     long add;
  632.     DSequence* aSeq= (DSequence*) Clone();
  633.     if (fromBases) add= StrLen( fromBases);  else add= 0;
  634.     if (add>0) {
  635.         char* h= aSeq->fBases;
  636.         char* m= aSeq->fMasks;
  637.         Boolean domask= (fMasks && fMasksOk);
  638.         long len = aSeq->fLength;
  639.         if (insertPoint<=0) {
  640.             fromBases= StrDup( fromBases);
  641.             StrExtendCat( &fromBases, h);
  642.             MemFree( h);
  643.             h= fromBases;
  644.             if (domask && fromMasks) {
  645.                 fromMasks= StrDup( fromMasks);
  646.                 StrExtendCat( &fromMasks, m);
  647.                 MemFree( m);
  648.                 m= fromMasks;
  649.                 }
  650.             }        
  651.         else if (insertPoint>=len) { 
  652.             StrExtendCat( &h, fromBases);
  653.             if (domask) StrExtendCat( &m, fromMasks);
  654.             }
  655.         else { 
  656.             h= (char*) MemMore( h, len+add);
  657.             MemMove( h+insertPoint+add, h+insertPoint, len-insertPoint);
  658.             MemMove( h+insertPoint+add, fromBases, add);
  659.             *((char*)h+len+add)= '\0';
  660.             if (domask) {
  661.                 m= (char*) MemMore( m, len+add);
  662.                 MemMove( m+insertPoint+add, m+insertPoint, len-insertPoint);
  663.                 MemMove( m+insertPoint+add, fromMasks, add);
  664.                 *((char*)m+len+add)= '\0';
  665.                 }
  666.             }
  667.         aSeq->fBases= h;
  668.         aSeq->fMasks= m;
  669.         }
  670.     aSeq->UpdateFlds();    
  671.     aSeq->Modified();
  672.     return aSeq;
  673. }
  674.  
  675.  
  676. DSequence* DSequence::Reverse()
  677. {
  678.     long i;
  679.     DSequence* aSeq= (DSequence*) Clone();    
  680.     if (fSelStart < 0) fSelStart= 0;
  681.     long len= aSeq->fLength - fSelStart;
  682.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  683.     char* pc= this->fBases+fSelStart + len;
  684.     char* pr= aSeq->fBases+fSelStart;
  685.     for (i= 0; i<len; i++) *pr++= *--pc;
  686.     Boolean domask= (fMasks && fMasksOk);
  687.     if (domask) {
  688.         pc= this->fMasks+fSelStart + len;
  689.         pr= aSeq->fMasks+fSelStart;
  690.         for (i= 0; i<len; i++) *pr++= *--pc;
  691.         }
  692.     aSeq->UpdateFlds();
  693.     aSeq->Modified();
  694.     return aSeq;
  695. }
  696.  
  697.  
  698. DSequence* DSequence::Complement()
  699. // A->T, T->A, G->C, C->G  
  700. {
  701.     Boolean        isrna;
  702.     DSequence* aSeq= (DSequence*) Clone();    
  703.     
  704.     if (fKind == kRNA) isrna= TRUE;
  705.     else if (fKind == kDNA || fKind == kNucleic) isrna= FALSE;
  706.     else return aSeq;    
  707.         
  708.     if (fSelStart < 0) fSelStart= 0;
  709.     long len= aSeq->fLength - fSelStart;
  710.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  711.     char* pc= this->fBases+fSelStart;
  712.     char* pr= aSeq->fBases+fSelStart;
  713.     NucleicComplement( isrna, pc, pr, len);
  714.         
  715.     aSeq->UpdateFlds();
  716.     aSeq->Modified();
  717.     return aSeq;
  718. }
  719.  
  720.  
  721.  
  722. DSequence* DSequence::ChangeCase(Boolean toUpper)
  723. // g->G or G->g
  724. {
  725.     DSequence* aSeq= (DSequence*) Clone();    
  726.     
  727.     if (fSelStart < 0) fSelStart= 0;
  728.     char* h= aSeq->fBases;
  729.     long len= aSeq->fLength - fSelStart;
  730.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  731.     char* pc= this->fBases+fSelStart;
  732.     char* pr= aSeq->fBases+fSelStart;
  733.     
  734.     long i;
  735.     if (toUpper) for (i= 0; i<len; i++) {
  736.         *pr++= toupper(*pc++);
  737.         }
  738.     else for (i= 0; i<len; i++) {
  739.         *pr++= tolower(*pc++);
  740.         }    
  741.         
  742.     aSeq->UpdateFlds();
  743.     aSeq->Modified();
  744.     return aSeq;
  745. }
  746.  
  747. DSequence* DSequence::Dna2Rna(Boolean toRna)
  748. // T->U or U->T 
  749. {
  750.     Boolean        isrna;
  751.     DSequence* aSeq= (DSequence*) Clone();    
  752.     
  753.     if (fKind==kRNA) isrna= TRUE;
  754.     else if (fKind == kDNA || fKind == kNucleic) isrna= FALSE;
  755.     else return aSeq;    
  756.      
  757.     if (fSelStart < 0) fSelStart= 0;
  758.     long len= aSeq->fLength - fSelStart;
  759.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  760.     char* pc= this->fBases+fSelStart;
  761.     char* pr= aSeq->fBases+fSelStart;
  762.     
  763.     long i;
  764.     char b, c;
  765.     if (toRna) for (i= 0; i<len; i++) {
  766.         b= *pc++;
  767.         if (b=='T') c= 'U';
  768.         else if (b=='t') c= 'u';
  769.         else c= b;
  770.         *pr++= c;
  771.         }
  772.     else for (i= 0; i<len; i++) {
  773.         b= *pc++;
  774.         if (b=='U') c= 'T';
  775.         else if (b=='u') c= 't';
  776.         else c= b;
  777.         *pr++= c;
  778.         }    
  779.         
  780.     aSeq->UpdateFlds();
  781.     aSeq->Modified();
  782.     return aSeq;
  783. }
  784.  
  785.  
  786. DSequence* DSequence::LockIndels( Boolean doLock)
  787. // indelSoft -> indelHard or vice versa 
  788. {    
  789.     char b, c, fromc, toc;
  790.     DSequence* aSeq= (DSequence*) Clone();    
  791.     
  792.     if (fSelStart < 0) fSelStart= 0;
  793.     long  len= aSeq->fLength - fSelStart;
  794.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  795.     char* pc= this->fBases+fSelStart;
  796.     char* pr= aSeq->fBases+fSelStart;
  797.     
  798.     if (doLock) { fromc= indelSoft; toc= indelHard; }
  799.     else { fromc= indelHard; toc= indelSoft; }
  800.     
  801.     for (long i= 0; i<len; i++) {
  802.         b= *pc++;
  803.         if (b == fromc) c= toc;
  804.         else c= b;
  805.         *pr++= c;
  806.         }    
  807.         
  808.     aSeq->UpdateFlds();
  809.     aSeq->Modified();
  810.     return aSeq;
  811. }
  812.  
  813.  
  814. DSequence* DSequence::Slide( long slideDist)
  815. /* insert "-" below (+dist) or above (-dist) of seq range,
  816.     squeeze out "-" above (below) of range
  817. */
  818. {     
  819.     long         slideSave, newlen, fulllen, len, i, j= 0;
  820.     char         b, cutchar, cutmask;
  821.     char        *pc, *pr;
  822.     
  823.     DSequence* aSeq= (DSequence*) Clone();    
  824.     if (fSelStart < 0) fSelStart= 0;
  825.     fulllen= aSeq->fLength;
  826.     len= fulllen - fSelStart;
  827.     Boolean domask= (fMasks && fMasksOk);
  828.     
  829.     if (slideDist < 0) {
  830.         // tough... do if below works...
  831.         // really need ability to extend seq in 5' (lower) direction 
  832.         //! build pr in reverse, then flip it....
  833.  
  834.         if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  835.         slideDist= -slideDist;
  836.         slideSave= slideDist;
  837.         aSeq->fBases= (char*) MemMore(aSeq->fBases, 2+fulllen+slideDist);
  838.         pr= aSeq->fBases;
  839.         pc= this->fBases + fulllen;
  840.         for (i= fulllen-1; i >= fSelStart+len; i--) *pr++= *--pc;
  841.              
  842.         if (pr == aSeq->fBases || *(pr-1) == indelEdge) cutchar= indelEdge;
  843.         else cutchar= indelSoft; //??
  844.         for (i=0; i<slideDist; i++) *pr++= cutchar;
  845.  
  846.         pc= this->fBases+fSelStart+len;
  847.         for (i= fSelStart+len-1; i >= fSelStart; i--) *pr++= *--pc;
  848.             
  849.         pc= this->fBases+fSelStart;
  850.         for (i= fSelStart-1; i>=0; i--) {
  851.             b= *--pc;
  852.             if (slideDist>0 && (b==indelSoft || b==indelEdge)) slideDist--;
  853.             else {
  854.                 *pr++= b;
  855.                 }
  856.             }
  857.  
  858.         newlen= pr - aSeq->fBases;
  859.         char* htmp= (char*) MemNew(newlen+1);
  860.         pc= htmp;
  861.         for (i=0; i<newlen; i++) *pc++= *--pr; //pr[j-i+1];
  862.         pc= '\0';
  863.         MemFree(aSeq->fBases);
  864.         aSeq->fBases = htmp;  
  865.         aSeq->fLength= newlen;
  866.         aSeq->fOrigin= fulllen - newlen; //== this->fLength - aSeq->fLength
  867.  
  868.         if (domask) {
  869.             slideDist= slideSave;
  870.             aSeq->fMasks= (char*) MemMore(aSeq->fMasks, 2+fulllen+slideDist);
  871.             pr= aSeq->fMasks;
  872.             pc= this->fMasks + fulllen;
  873.             for (i= fulllen-1; i >= fSelStart+len; i--) *pr++= *--pc;
  874.                  
  875.             cutmask= 0x80; // required empty mask  
  876.             for (i=0; i<slideDist; i++) *pr++= cutmask;
  877.     
  878.             pc= this->fMasks+fSelStart+len;
  879. #ifdef WIN_MSWIN
  880.             // mswin is crashing here for no obvious reason
  881.             //pc--; len--;
  882.             for (i= fSelStart + len - 1; i >= fSelStart; i--) {
  883.                  register short amask;
  884.                  --pc;
  885.                  amask= *pc;    //mswin bomb is here !
  886.                  *pr= amask;
  887.                  pr++;
  888.          }
  889. #else
  890.             for (i= fSelStart+len-1; i >= fSelStart; i--) *pr++= *--pc;
  891. #endif
  892.             pc= this->fMasks+fSelStart;
  893.             for (i= fSelStart-1; i>=0; i--) {
  894.                 b= *--pc;
  895.                 if (slideDist>0 && (b==cutmask)) slideDist--;
  896.                 else {
  897.                     *pr++= b;
  898.                     }
  899.                 }
  900.                 
  901.             newlen= pr - aSeq->fMasks;
  902.             htmp= (char*) MemNew(newlen+1);
  903.             pc= htmp;
  904.             for (i=0; i<newlen; i++) *pc++= *--pr; //pr[j-i+1];
  905.             pc= '\0';
  906.             MemFree(aSeq->fMasks);
  907.             aSeq->fMasks = htmp;  
  908.             }
  909.             
  910.         }    
  911.         
  912.     else if (slideDist > 0) {
  913.         fulllen= len;
  914.         if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  915.         aSeq->fBases= (char*) MemMore(aSeq->fBases, 2+fSelStart+fulllen+slideDist); //make it large enough
  916.  
  917.         pr= aSeq->fBases;
  918.       if (fSelStart==0 || pr[fSelStart-1]==indelEdge) cutchar= indelEdge; 
  919.         else cutchar= indelSoft; 
  920.  
  921.         pc= this->fBases+fSelStart;
  922.         pr= aSeq->fBases+fSelStart;
  923.  
  924.         for (j= 0; j<slideDist; j++) *pr++= cutchar;
  925.         for (i= 0; i<len; i++) *pr++= *pc++; 
  926.     
  927.         for (i= len; i<fulllen; i++) {
  928.             b= *pc++;
  929.             if (slideDist>0 && (b==indelSoft || b==indelEdge)) 
  930.                 slideDist--;
  931.             else {
  932.                 *pr++= b;
  933.                 }
  934.             }
  935.         
  936.         *pr= 0;
  937.         newlen= pr - aSeq->fBases; 
  938.         aSeq->fBases= (char*)MemMore(aSeq->fBases, newlen+1);  
  939.         
  940.         if (domask) {
  941.             aSeq->fMasks= (char*) MemMore(aSeq->fMasks, 2+fSelStart+fulllen+slideDist);
  942.           cutmask= 0x80; 
  943.             pc= this->fMasks+fSelStart;
  944.             pr= aSeq->fMasks+fSelStart;
  945.     
  946.             for (j= 0; j<slideDist; j++) *pr++= cutmask;
  947.             for (i= 0; i<len; i++) *pr++= *pc++; 
  948.  
  949.             for (i= len; i<fulllen; i++) {
  950.                 b= *pc++;
  951.                 if (slideDist>0 && (b==cutmask)) // this is BAD? must be done w/ base slide !?
  952.                     slideDist--;
  953.                 else {
  954.                     *pr++= b;
  955.                     }
  956.                 }
  957.             
  958.             *pr= 0;
  959.             newlen= pr - aSeq->fMasks; 
  960.             aSeq->fMasks= (char*)MemMore(aSeq->fMasks, newlen+1);  
  961.             }
  962.         }
  963.             
  964.     aSeq->UpdateFlds();
  965.     aSeq->Modified();
  966.     return aSeq;
  967. }
  968.  
  969.  
  970. DSequence* DSequence::Compress()
  971. // pull out "-", ".", ? anything but nuc/aa codes ? 
  972. {
  973.     long  len, fulllen, i;
  974.     char     b, c;
  975.     char * mc, * mr;
  976.     Boolean domask= (fMasks && fMasksOk);
  977.  
  978.     DSequence* aSeq= (DSequence*) Clone();    
  979.     if (fSelStart < 0) fSelStart= 0;
  980.     len= aSeq->fLength - fSelStart;
  981.     fulllen= len;
  982.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  983.     char* pc= this->fBases+fSelStart;
  984.     char* pr= aSeq->fBases+fSelStart;
  985.     if (domask) {
  986.         mc= this->fMasks+fSelStart;
  987.         mr= aSeq->fMasks+fSelStart;
  988.         }
  989.         
  990.     if (domask) for (i= 0; i<len; i++) {
  991.         b= *pc++;
  992.         c= *mc++;
  993.         if (b!=indelSoft && b!=indelEdge) { 
  994.             *pr++= b;
  995.             *mr++= c;
  996.             }
  997.         }
  998.     else for (i= 0; i<len; i++) {
  999.         b= *pc++;
  1000.         if (b!=indelSoft && b!=indelEdge) 
  1001.             *pr++= b;
  1002.         }
  1003.     for (i=len; i<fulllen; i++) *pr++= *pc++;
  1004.     if (domask) for (i=len; i<fulllen; i++) *mr++= *mc++;
  1005.  
  1006.     long newlen= pr - aSeq->fBases;
  1007.     aSeq->fBases= (char*) MemMore( aSeq->fBases, newlen+1);
  1008.     aSeq->fBases[newlen]= 0;
  1009.     if (domask) {
  1010.         aSeq->fMasks= (char*) MemMore( aSeq->fMasks, newlen+1);
  1011.         aSeq->fMasks[newlen]= 0;
  1012.         }
  1013.         
  1014.     aSeq->UpdateFlds();
  1015.     aSeq->Modified();
  1016.     return aSeq;
  1017. }
  1018.  
  1019.  
  1020. DSequence* DSequence::Translate(Boolean keepUnselected)
  1021. {    
  1022.     char  c;
  1023.     char    acodon[3];
  1024.     Boolean isrna, nocodon;
  1025.     char    *pr, *pc, *mr, *mc;
  1026.     long    len, fulllen, newlen, i, j, k;
  1027.     
  1028.     if (DCodons::NotAvailable()) return NULL;
  1029.     DSequence* aSeq= (DSequence*) Clone();        
  1030.     // !? what do we do w/ fMasks !?
  1031.     Boolean domask= (fMasks && fMasksOk);
  1032.     
  1033.     if (fSelStart < 0) fSelStart= 0;
  1034.     fulllen= this->fLength - fSelStart;
  1035.     len= fulllen;
  1036.     if (fSelBases > 0 && fSelBases < len) len= fSelBases;
  1037.     
  1038.     switch (fKind) {
  1039.     
  1040.         case kRNA:
  1041.         case kDNA:
  1042.         case kNucleic: 
  1043.             {
  1044.             if (keepUnselected) {
  1045.                 pr= aSeq->fBases + fSelStart; // start at selection, leaving untranslated 5'
  1046.                 if (domask) mr= aSeq->fMasks + fSelStart;
  1047.                 }
  1048.             else {
  1049.                 pr= aSeq->fBases;
  1050.                  if (domask) mr= mr= aSeq->fMasks;
  1051.                  }
  1052.                  
  1053.             pc= this->fBases + fSelStart;
  1054.             if (domask) mc= this->fMasks + fSelStart;
  1055.             isrna= (fKind == kRNA);
  1056.             i= fSelStart;
  1057.             while (i < len) {
  1058.                 for (j=0; j<3; j++) {
  1059.                     c= toupper( *pc++);
  1060.                     if (isrna && c=='U') c= 'T';
  1061.                     acodon[j]= c;
  1062.                     }
  1063.                 i += 3;
  1064.                 for (k= 0, nocodon= true; k<64; k++) 
  1065.                     if ( MemCmp(DCodons::codon(k),acodon,3) == 0) { 
  1066.                         *pr++= DCodons::amino(k); //gCodonTable[k].amino;
  1067.                         nocodon= false;
  1068.                         break;
  1069.                         }
  1070.                 if (nocodon) *pr++= '?';
  1071.                 if (domask) {
  1072.                     *mr++= *mc;
  1073.                     mc += 3;
  1074.                     }
  1075.                 }
  1076.                  // copy over untranslated portion
  1077.             if (keepUnselected) {
  1078.                 long savei= i;
  1079.                 for ( ; i < fulllen; i++)  *pr++= *pc++; 
  1080.                 if (domask) for ( i= savei; i < fulllen; i++)  *mr++= *mc++; 
  1081.                 }
  1082.             }
  1083.             break;
  1084.             
  1085.         case kAmino: 
  1086.             {
  1087.             char *bestcodon;
  1088.             long  selend;
  1089.             
  1090.             if (len < fulllen) {
  1091.                 newlen= fulllen + len*3;
  1092.                 selend= fSelStart + len;
  1093.                 }
  1094.             else {
  1095.                 newlen= fulllen * 3;
  1096.                 selend= len;
  1097.                 }
  1098.             pc= this->fBases + fSelStart;
  1099.             aSeq->fBases= (char*) MemMore( aSeq->fBases, newlen);
  1100.             if (keepUnselected)
  1101.                 pr= aSeq->fBases + fSelStart;
  1102.             else
  1103.                 pr= aSeq->fBases;
  1104.                 
  1105.             if (domask) {
  1106.                 mc= this->fMasks + fSelStart;
  1107.                 aSeq->fMasks= (char*) MemMore( aSeq->fMasks, newlen);
  1108.                 if (keepUnselected)
  1109.                     mr= aSeq->fMasks + fSelStart;
  1110.                 else
  1111.                     mr= aSeq->fMasks;
  1112.                 }
  1113.                 
  1114.             for (i= fSelStart; i<selend; i++) {
  1115.                 c= toupper( *pc++);
  1116.                 bestcodon= (char*) DCodons::FindBestCodon(c);
  1117.                 MemCpy( pr, bestcodon, 3);  
  1118.                 pr += 3;
  1119.                 if (domask) {
  1120.                     for (j=0; j<3; j++) *mr++= *mc;
  1121.                     mc++;
  1122.                     }
  1123.                 }
  1124.                  // copy over untranslated portion
  1125.             if (keepUnselected) {
  1126.                 long savei= i;
  1127.                 for ( ; i < fulllen; i++)  *pr++= *pc++; 
  1128.                 if (domask) for ( i=savei; i < fulllen; i++)  *mr++= *mc++; 
  1129.                 }
  1130.             }
  1131.             break;
  1132.         
  1133.         default:
  1134.             return NULL;
  1135.          
  1136.       }
  1137.         
  1138.     *pr= '\0';
  1139.     newlen= pr - aSeq->fBases;
  1140.     aSeq->fBases= (char*) MemMore( aSeq->fBases, newlen);
  1141.     if (domask) {
  1142.         *mr= '\0';
  1143.         newlen= mr - aSeq->fMasks;
  1144.         aSeq->fMasks= (char*) MemMore( aSeq->fMasks, newlen);
  1145.         }
  1146.     aSeq->UpdateFlds();
  1147.     aSeq->Modified();
  1148.     return aSeq;    
  1149. }
  1150.  
  1151.  
  1152. //Boolean domask= (fMasks && fMasksOk);
  1153.  
  1154.  
  1155. void DSequence::DoWrite(DFile* aFile, short format)
  1156. {  
  1157.     if (fLength>0) {
  1158.         DSeqFile::WriteSeqWrapper( aFile, fBases, fLength, format, fName);
  1159.         if (fMasks && fMasksOk && DSeqFile::gWriteMasks) 
  1160.             DSeqFile::WriteMaskWrapper( aFile, fMasks, fLength, format, fName);
  1161.         }
  1162. }
  1163.  
  1164.  
  1165. void DSequence::DoWriteSelection(DFile* aFile, short format)
  1166. {
  1167.     if (fSelBases==0) 
  1168.         DoWrite( aFile, format);
  1169.     else if (fLength>0) {
  1170.         long aStart= Min( fSelStart, fLength);
  1171.         long aLength= Min( fSelBases, fLength - fSelStart);
  1172.         DSeqFile::WriteSeqWrapper(  aFile, fBases+aStart, aLength, format, fName);
  1173.         if (fMasks && fMasksOk && DSeqFile::gWriteMasks)  
  1174.             DSeqFile::WriteMaskWrapper( aFile, fMasks+aStart, aLength, format, fName);
  1175.         }
  1176. }
  1177.  
  1178.  
  1179.  
  1180.  
  1181. Boolean DSequence::GoodChar(char& aChar) 
  1182. {
  1183.     char* aacodes = gAminos;
  1184.     char* iupaccodes = gIubbase;
  1185.     char* nacodes = gPrimenuc;
  1186.     char    c;
  1187.     char* codes = NULL;
  1188.  
  1189.   if (aChar=='\n' || aChar=='\r')  
  1190.         return FALSE;
  1191.     else if (aChar<' ' || aChar>'~')  
  1192.         return TRUE; //allow return,tab,delete, ? option chars
  1193.  
  1194.   switch (fKind) {
  1195.       case kDNA            :
  1196.         case kRNA            : codes= nacodes; break;
  1197.         case kNucleic    : codes= iupaccodes; break;
  1198.         case kAmino        : codes= aacodes; break;
  1199.         default:
  1200.       case kOtherSeq: return TRUE;
  1201.         }
  1202.         
  1203.   c= toupper(aChar);
  1204.     while (*codes) {
  1205.         if (c == *codes) {
  1206.           switch (fKind) {
  1207.                 case kDNA: if (c == 'U') {
  1208.                     if (aChar=='u') aChar= 't'; else aChar= 'T';
  1209.                     }
  1210.                 break;
  1211.                 case kRNA: if (c == 'T') {
  1212.                     if (aChar=='t') aChar= 'u'; else aChar= 'U';
  1213.                     }
  1214.                 break;
  1215.                 }
  1216.             return TRUE;
  1217.             }
  1218.         codes++;
  1219.         }
  1220.     return FALSE;
  1221. }
  1222.  
  1223.  
  1224. Boolean DSequence::IsConsensus()
  1225. {
  1226.     return (fName && StringCmp(fName,kConsensus)==0);
  1227. }
  1228.  
  1229.  
  1230.  
  1231. DSequence* MakeSequence(char* name, char* bases, char* info, long modtime)
  1232. {
  1233.     DSequence* aSeq= new DSequence();
  1234.     aSeq->SetName(name);
  1235.     //?? can we always assume bases is null terminated?
  1236.     aSeq->SetBases(bases);
  1237.     aSeq->SetInfo(info);
  1238.      if (modtime!=0) aSeq->SetTime( (unsigned long)modtime);
  1239.     aSeq->UpdateFlds();      
  1240.     return aSeq;
  1241. }
  1242.  
  1243.  
  1244.  
  1245.  
  1246.  
  1247. long DSequence::SearchAgain() 
  1248. {
  1249.     return Search( "", FALSE);
  1250. }
  1251.  
  1252.  
  1253.  
  1254. void DSequence::SearchORF( long& start, long& stop) 
  1255. {
  1256.     char  startc[2];
  1257.     long  i, j, k;
  1258.     short    ns;
  1259.     CodonStat  stops[10];
  1260.     
  1261.     start= kSearchNotFound;
  1262.     stop = kSearchNotFound;
  1263.     if (DCodons::NotAvailable()) return;
  1264.     
  1265.     for (ns= 0, i= 0; i<64; i++) {
  1266.         if (ns<10 && DCodons::amino(i) == '*') 
  1267.             stops[ns++] = DCodons::fCodons[i];
  1268.         }
  1269.     if (ns == 0) return;
  1270.         
  1271.     switch (fKind) {  
  1272.         case kRNA:
  1273.         case kDNA:
  1274.         case kNucleic: 
  1275.              start= Search( DCodons::startcodon(), FALSE);
  1276.              if (start == kSearchNotFound) return;
  1277.              for (j= start+3; j<fLength; j += 3) {
  1278.                  for (k=0; k<ns; k++) {
  1279.                      if (StrNICmp( fBases+j, stops[k].codon, 3) == 0) {
  1280.                          stop= j+3; // +3 to include this codon
  1281.                          return;
  1282.                          }
  1283.                      }
  1284.                  }
  1285.             break;
  1286.             
  1287.         case kAmino:
  1288.             startc[0]= DCodons::startamino();
  1289.             startc[1]= 0;
  1290.              start= Search( startc, FALSE);
  1291.              if (start == kSearchNotFound) return;
  1292.              for (j= start+1; j<fLength; j++) {
  1293.                  for (k=0; k<ns; k++) {
  1294.                      if (toupper(fBases[j]) == toupper(stops[k].amino)) {
  1295.                          stop= j+1; // +1 to include this codon
  1296.                          return;
  1297.                          }
  1298.                      }
  1299.                  }
  1300.             break;
  1301.             
  1302.         default:
  1303.             break;
  1304.         }
  1305. }
  1306.  
  1307.  
  1308.  
  1309. void DSequence::NucleicSearch(char* source, char* target, SearchRec& aSearchRec)
  1310. {
  1311.     long k, j, tBits, lens, ti;
  1312.     Boolean done;
  1313.  
  1314.                     //find target[ti] where tBits == kMaskA,C,G,T for fastest search
  1315.                     //This part is wasted time for SearchAgain... store tBits/ti in saverec...
  1316.     lens= StrLen(target);
  1317.     for (j= 0, ti= 0, tBits= 0; j<lens && !tBits; j++) {
  1318.         ti= j;
  1319.         tBits= NucleicBits( target[ti]);
  1320.         tBits &= kMaskNucs; //    bClr( tBits, kBitExtra); //drop RNA bit
  1321.         }
  1322.  
  1323.     do {
  1324. reloop:
  1325.             aSearchRec.indx += aSearchRec.step; 
  1326.             aSearchRec.lim  -= aSearchRec.step;    
  1327.             
  1328.             if (aSearchRec.step>0) {
  1329.                 for (j= 0; j<aSearchRec.lim; j++)
  1330.                  // if (IsNucleicMatch( tBits, NucleicBits( source[j+aSearchRec.indx]))) 
  1331.                  if ( tBits & NucleicBits( source[j+aSearchRec.indx]) ) 
  1332.                  {
  1333.                     k= j; 
  1334.                     goto nextindx;
  1335.                  }
  1336.                 k= aSearchRec.lim;
  1337.                 }            
  1338.             else {  //step<0 IS BAD .. see UTextDoc version...
  1339.                 for (j= 0; j > -aSearchRec.lim; j--) 
  1340.                   //if (IsNucleicMatch( tBits, NucleicBits( source[j+aSearchRec.indx])))
  1341.                   if ( tBits & NucleicBits( source[j+aSearchRec.indx]) ) 
  1342.                 {
  1343.                     k= j; 
  1344.                     goto nextindx;
  1345.                     }
  1346.                 k= -aSearchRec.lim;
  1347.                 }
  1348.                 
  1349. nextindx:  
  1350.             aSearchRec.indx += k; 
  1351.             aSearchRec.lim  -= k;
  1352.                 //? not enough source left to match full target 
  1353.             if (aSearchRec.step>0) {
  1354.                 if (aSearchRec.lim < lens-ti) aSearchRec.lim= 0;
  1355.                 }            
  1356.             else {
  1357.                 if (aSearchRec.lim < ti) aSearchRec.lim= 0; //? or lim < ti-1 
  1358.                 }
  1359.             done= (aSearchRec.lim==0);
  1360.             
  1361.             if (!done) {
  1362.                 long srcstart= aSearchRec.indx - ti;
  1363.                 for (j= 0; j<lens; j++) 
  1364.                  if (j != ti) {
  1365.                     //if (!IsNucleicMatch( NucleicBits( target[j]), NucleicBits( source[j+srcstart]))) 
  1366.                   if (!(kMaskNucs & NucleicBits( target[j]) & NucleicBits( source[j+srcstart]))) 
  1367.                       goto reloop;
  1368.                     }
  1369.                 done= true;
  1370.                 }
  1371.         } while (!done);
  1372. }  
  1373.     
  1374.  
  1375. void DSequence::ProteinSearch(char* source, char* target, SearchRec& aSearchRec)
  1376. {
  1377.     //just toupper match of chars 
  1378.     char      tChar; 
  1379.     Boolean  done;
  1380.     long        lens, k, j;
  1381.  
  1382.     lens= StrLen( target);
  1383.     tChar= toupper( target[0]);
  1384.     
  1385.     do {
  1386. reloop:
  1387.             aSearchRec.indx += aSearchRec.step; 
  1388.             aSearchRec.lim  -= aSearchRec.step;    
  1389.             
  1390.             if (aSearchRec.step>0) {
  1391.                 for (j= 0; j<aSearchRec.lim; j++) 
  1392.                  if (tChar == toupper( source[j+aSearchRec.indx])) {
  1393.                     k= j;
  1394.                     goto nextindx;
  1395.                     }
  1396.                 k= aSearchRec.lim;
  1397.                 }            
  1398.             else {
  1399.                 for (j= 0; j > -aSearchRec.lim; j--)
  1400.                  if (tChar == toupper(source[j+aSearchRec.indx])) {
  1401.                     k= j; 
  1402.                     goto nextindx;
  1403.                     }
  1404.                 k= -aSearchRec.lim;
  1405.                 }
  1406.                 
  1407. nextindx:  
  1408.             aSearchRec.indx += k; 
  1409.             aSearchRec.lim  -= k;
  1410.                 //? not enough source left to match full target 
  1411.             if (aSearchRec.step>0) {
  1412.                 if (aSearchRec.lim < lens-1) aSearchRec.lim= 0;
  1413.                 }            
  1414.             else {
  1415.                 if (aSearchRec.lim < 1) aSearchRec.lim= 0;  
  1416.                 }
  1417.             done= (aSearchRec.lim==0);
  1418.             
  1419.             if (!done) {
  1420.                 for (j= 1; j<lens; j++) {
  1421.                     if (toupper( target[j]) != toupper( source[j+aSearchRec.indx])) 
  1422.                         goto reloop;
  1423.                     }
  1424.                 done= true;
  1425.                 }
  1426.         } while (!done);
  1427. }  
  1428.  
  1429.  
  1430.     
  1431.  
  1432. long DSequence::Search( char* target, Boolean backwards)
  1433. /* search for target in sequence;
  1434.     search == abs index into fBases[0..fLength] or kSearchNotFound
  1435.   if (target='') then find next  
  1436.   if (backwards) then backwards search
  1437. */
  1438. {
  1439.     long        limit, aStart;
  1440.     SearchRec  aSearchRec;
  1441.  
  1442.     aStart= fSelStart;
  1443.     if (!target || !*target) {  //search again
  1444.         if (fSearchRec.targ && (fSearchRec.step == -1 || fSearchRec.step == 1)) {  
  1445.             aSearchRec= fSearchRec;
  1446.             target= aSearchRec.targ;
  1447.             }        
  1448.         else {
  1449.             return kSearchNotFound;
  1450.             }
  1451.         }        
  1452.     else {
  1453.         limit= fLength; //StrLen(fBases);
  1454.         if (aStart < 0) aStart= 0;
  1455.         if (backwards) {
  1456.             if (!aStart) aStart= limit-1;
  1457.             else limit= aStart+1;  //? off-by-one ??
  1458.             }        
  1459.         else {
  1460.             if (aStart > limit) aStart= 0; 
  1461.             else limit -= aStart;
  1462.             }
  1463.         if (fSelBases > 0 && fSelBases < limit) limit= fSelBases;
  1464.         if (backwards) aSearchRec.step= -1;
  1465.         else aSearchRec.step= 1;
  1466.         aSearchRec.indx= aStart - aSearchRec.step;  //indx always in [0..handsize-1], except 1st is -1 
  1467.         aSearchRec.lim = limit + aSearchRec.step;  //lim always in [1..handsize] or 0=done, except 1st is hand+1 
  1468.         aSearchRec.targ= StrDup(target);
  1469.         }
  1470.  
  1471.     switch (fKind) {  
  1472.         case kRNA:
  1473.         case kDNA:
  1474.         case kNucleic: 
  1475.             NucleicSearch(fBases,target,aSearchRec); 
  1476.             break;
  1477.         default:
  1478.             ProteinSearch(fBases,target,aSearchRec);
  1479.             break;
  1480.         }
  1481.     
  1482.     fSearchRec= aSearchRec;
  1483.     if (aSearchRec.lim == 0) return kSearchNotFound;
  1484.     else return aSearchRec.indx;
  1485. }  
  1486.  
  1487.  
  1488.